library(swimplot) library(grid) library(gtable) library(readr) library(mosaic) library(dplyr) library(survival) library(survminer) library(ggplot2) library(scales) library(coxphf) library(ggthemes) library(tidyverse) library(gtsummary) library(flextable) library(parameters) library(car) library(ComplexHeatmap) library(tidyverse) library(readxl) library(janitor) library(DT) library(pROC) library(rms)
#ctDNA Detection Rates by Window and Stages
#ctDNA at Baseline
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Base!="",]
circ_data$ctDNA.Base <- factor(circ_data$ctDNA.Base, levels=c("NEGATIVE","POSITIVE"))
circ_data <- subset(circ_data, ctDNA.Base %in% c("NEGATIVE", "POSITIVE"))
circ_data$Stage <- factor(circ_data$Stage, levels=c("O","I","II","III","IV"))
positive_counts_by_stage <- aggregate(circ_data$ctDNA.Base == "POSITIVE", by=list(circ_data$Stage), FUN=sum)
total_counts_by_stage <- aggregate(circ_data$ctDNA.Base, by=list(circ_data$Stage), FUN=length)
combined_data <- data.frame(
Stage = total_counts_by_stage$Group.1,
Total_Count = total_counts_by_stage$x,
Positive_Count = positive_counts_by_stage$x,
Rate = (positive_counts_by_stage$x / total_counts_by_stage$x) * 100 # Convert to percentage
)
combined_data$Rate <- sprintf("%.2f%%", combined_data$Rate)
overall_total_count <- nrow(circ_data)
overall_positive_count <- nrow(circ_data[circ_data$ctDNA.Base == "POSITIVE",])
overall_positivity_rate <- (overall_positive_count / overall_total_count) * 100 # Convert to percentage
overall_row <- data.frame(
Stage = "Overall",
Total_Count = overall_total_count,
Positive_Count = overall_positive_count,
Rate = sprintf("%.2f%%", overall_positivity_rate)
)
combined_data <- rbind(combined_data, overall_row)
print(combined_data)
#ctDNA post-treatment
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.postTx!="",]
circ_data$ctDNA.postTx <- factor(circ_data$ctDNA.postTx, levels=c("NEGATIVE","POSITIVE"))
circ_data$Stage <- factor(circ_data$Stage, levels=c("O","I","II","III","IV"))
positive_counts_by_stage <- aggregate(circ_data$ctDNA.postTx == "POSITIVE", by=list(circ_data$Stage), FUN=sum)
total_counts_by_stage <- aggregate(circ_data$ctDNA.postTx, by=list(circ_data$Stage), FUN=length)
combined_data <- data.frame(
Stage = total_counts_by_stage$Group.1,
Total_Count = total_counts_by_stage$x,
Positive_Count = positive_counts_by_stage$x,
Rate = (positive_counts_by_stage$x / total_counts_by_stage$x) * 100 # Convert to percentage
)
combined_data$Rate <- sprintf("%.2f%%", combined_data$Rate)
overall_total_count <- nrow(circ_data)
overall_positive_count <- nrow(circ_data[circ_data$ctDNA.postTx == "POSITIVE",])
overall_positivity_rate <- (overall_positive_count / overall_total_count) * 100 # Convert to percentage
overall_row <- data.frame(
Stage = "Overall",
Total_Count = overall_total_count,
Positive_Count = overall_positive_count,
Rate = sprintf("%.2f%%", overall_positivity_rate)
)
combined_data <- rbind(combined_data, overall_row)
print(combined_data)
#Overview plot
setwd("~/Downloads")
clinstage <- read.csv("EORTC ICI_OP.csv")
clinstage_df <- as.data.frame(clinstage)
# Creating the basic swimmer plot
oplot <- swimmer_plot(df=clinstage_df,
id='PatientName',
end='fu.diff.months',
fill='gray',
width=.01)
# Adding themes and scales
oplot <- oplot + theme(panel.border = element_blank())
oplot <- oplot + scale_y_continuous(breaks = seq(0, 48, by = 3))
oplot <- oplot + labs(x ="Patients", y="Months from Immunotherapy Start")
# Adding swimmer points
oplot_ev1 <- oplot + swimmer_points(df_points=clinstage_df,
id='PatientName',
time='date.diff.months',
name_shape ='Event_type',
name_col = 'Event',
size=3.5,fill='black')
# Optionally uncomment and use col='darkgreen' if needed
# Adding shape manual scale
oplot_ev1.1 <- oplot_ev1 + ggplot2::scale_shape_manual(name="Event_type",
values=c(1,16,6,4),
breaks=c('ctDNA_neg','ctDNA_pos','Imaging','Death'))
# Display the plot
oplot_ev1.1
oplot_ev2 <- oplot_ev1.1 + swimmer_lines(df_lines=clinstage_df,
id='PatientName',
start='Tx_start.months',
end='Tx_end.months',
name_col='Tx_type',
size=3.5,
name_alpha = 1.0)
oplot_ev2 <- oplot_ev2 + guides(linetype = guide_legend(override.aes = list(size = 5, color = "black")))
oplot_ev2
oplot_ev2.2 <- oplot_ev2 + ggplot2::scale_color_manual(name="Event",values=c( "orange", "black", "black", "lightblue", "green", "red"))
oplot_ev2.2
#Overview plot - stratified by BOR
setwd("~/Downloads")
clinstage <- read.csv("EORTC ICI_OP.csv")
clinstage_df <- as.data.frame(clinstage)
# Creating the basic swimmer plot
oplot_stratify <-swimmer_plot(df=clinstage_df,
id='PatientName',
end='fu.diff.months',
col="gray",
alpha=0.75,
width=.01,
base_size = 14,
stratify= c('RECIST'))
oplot_stratify <- oplot_stratify + theme(panel.border = element_blank())
oplot_stratify <- oplot_stratify + scale_y_continuous(breaks = seq(0, 42, by = 3))
oplot_stratify <- oplot_stratify + labs(x ="Patients" , y="Months from Immunotherapy Start")
# Adding swimmer points
oplot_ev1 <- oplot_stratify + swimmer_points(df_points=clinstage_df,
id='PatientName',
time='date.diff.months',
name_shape ='Event_type',
name_col = 'Event',
size=3.5,fill='black')
# Optionally uncomment and use col='darkgreen' if needed
# Adding shape manual scale
oplot_ev1.1 <- oplot_ev1 + ggplot2::scale_shape_manual(name="Event_type",
values=c(1,16,6,4),
breaks=c('ctDNA_neg','ctDNA_pos','Imaging','Death'))
# Display the plot
oplot_ev1.1
oplot_ev2 <- oplot_ev1.1 + swimmer_lines(df_lines=clinstage_df,
id='PatientName',
start='Tx_start.months',
end='Tx_end.months',
name_col='Tx_type',
size=3.5,
name_alpha = 1.0)
oplot_ev2 <- oplot_ev2 + guides(linetype = guide_legend(override.aes = list(size = 5, color = "black")))
oplot_ev2
oplot_ev2.2 <- oplot_ev2 + ggplot2::scale_color_manual(name="Event",values=c( "orange", "black", "black", "lightblue", "green", "red"))
oplot_ev2.2
#Association of Baseline ctDNA MTM levels with clinicopathological factors
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Base!="",]
circ_datadf <- as.data.frame(circ_data)
tally(~cStage, data=circ_data, margins = TRUE)
cStage
0-II III-IV Total
7 22 29
circ_data$cStage <- factor(circ_data$cStage, levels = c("0-II","III-IV"), labels = c("0-II (n=7)","III-IV (n=22)"))
boxplot(ctDNA.Base.MTM~cStage, data=circ_data, main="ctDNA pre-treatment MTM - Stage", xlab="Stage", ylab="MTM/mL", col="white",border="black")
m1<-wilcox.test(ctDNA.Base.MTM ~ cStage, data=circ_data, na.rm=TRUE, exact=FALSE, conf.int=TRUE)
print(m1)
Wilcoxon rank sum test with continuity correction
data: ctDNA.Base.MTM by cStage
W = 39.5, p-value = 0.05896
alternative hypothesis: true location shift is not equal to 0
95 percent confidence interval:
-5.19000e+01 4.27605e-05
sample estimates:
difference in location
-6.799975
tally(~Inclusion.status, data=circ_data, margins = TRUE)
Inclusion.status
Loco-regional Metastatic Total
18 11 29
circ_data$Inclusion.status <- factor(circ_data$Inclusion.status, levels = c("Loco-regional","Metastatic"), labels = c("Loco-regional (n=18)","Metastatic (n=11)"))
boxplot(ctDNA.Base.MTM~Inclusion.status, data=circ_data, main="ctDNA pre-treatment MTM - Disease Status", xlab="Disease Status", ylab="MTM/mL", col="white",border="black")
m2<-wilcox.test(ctDNA.Base.MTM ~ Inclusion.status, data=circ_data, na.rm=TRUE, exact=FALSE, conf.int=TRUE)
print(m2)
Wilcoxon rank sum test with continuity correction
data: ctDNA.Base.MTM by Inclusion.status
W = 93.5, p-value = 0.8219
alternative hypothesis: true location shift is not equal to 0
95 percent confidence interval:
-8.399949 6.899953
sample estimates:
difference in location
-0.3727674
tally(~cT.Status, data=circ_data, margins = TRUE)
cT.Status
cT0-T2 cT3-T4 cTx Total
11 17 1 29
circ_data$cT.Status <- factor(circ_data$cT.Status, levels = c("cT0-T2","cT3-T4"), labels = c("cT0-T2 (n=11)","cT3-T4 (n=17)"))
boxplot(ctDNA.Base.MTM~cT.Status, data=circ_data, main="ctDNA pre-treatment MTM - cT status", xlab="cT status", ylab="MTM/mL", col="white",border="black")
m3<-wilcox.test(ctDNA.Base.MTM ~ cT.Status, data=circ_data, na.rm=TRUE, exact=FALSE, conf.int=TRUE)
print(m3)
Wilcoxon rank sum test with continuity correction
data: ctDNA.Base.MTM by cT.Status
W = 45.5, p-value = 0.02523
alternative hypothesis: true location shift is not equal to 0
95 percent confidence interval:
-57.0000356 -0.7999756
sample estimates:
difference in location
-7.400021
tally(~cN.Status, data=circ_data, margins = TRUE)
cN.Status
cN0 cN1-N3 cNx Total
12 15 2 29
circ_data$cN.Status <- factor(circ_data$cN.Status, levels = c("cN0","cN1-N3"), labels = c("cN0 (n=12)","cN1-N3 (n=15)"))
boxplot(ctDNA.Base.MTM~cN.Status, data=circ_data, main="ctDNA pre-treatment MTM - cN status", xlab="cN status", ylab="MTM/mL", col="white",border="black")
m4<-wilcox.test(ctDNA.Base.MTM ~ cN.Status, data=circ_data, na.rm=TRUE, exact=FALSE, conf.int=TRUE)
print(m4)
Wilcoxon rank sum test with continuity correction
data: ctDNA.Base.MTM by cN.Status
W = 32.5, p-value = 0.005379
alternative hypothesis: true location shift is not equal to 0
95 percent confidence interval:
-60.200006 -3.299997
sample estimates:
difference in location
-10.15193
tally(~p16.status, data=circ_data, margins = TRUE)
p16.status
Negative Positive Unknown Total
24 4 1 29
circ_data$p16.status <- factor(circ_data$p16.status, levels = c("Negative","Positive"), labels = c("p16 neg (n=24)","p16 pos (n=4)"))
boxplot(ctDNA.Base.MTM~p16.status, data=circ_data, main="ctDNA pre-treatment MTM - p16 status", xlab="p16 status", ylab="MTM/mL", col="white",border="black")
m5<-wilcox.test(ctDNA.Base.MTM ~ p16.status, data=circ_data, na.rm=TRUE, exact=FALSE, conf.int=TRUE)
print(m5)
Wilcoxon rank sum test with continuity correction
data: ctDNA.Base.MTM by p16.status
W = 55.5, p-value = 0.6453
alternative hypothesis: true location shift is not equal to 0
95 percent confidence interval:
-47.89994 46.29998
sample estimates:
difference in location
0.8901786
#Median MTM/mL levels for ctDNA positive pts pre-treatment
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Base!="",]
circ_data <- circ_data[circ_data$ctDNA.Base=="POSITIVE",]
median_ctDNA <- median(circ_data$ctDNA.Base.MTM, na.rm = TRUE)
range_ctDNA <- range(circ_data$ctDNA.Base.MTM, na.rm = TRUE)
cat("Median MTM/mL post-treatment:", median_ctDNA, "\n")
Median MTM/mL post-treatment: 8.4
cat("Range MTM/mL post-treatment:", range_ctDNA[1], "-", range_ctDNA[2], "\n")
Range MTM/mL post-treatment: 0.2 - 986
#Median MTM/mL levels for ctDNA positive pts post-treatment
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.postTx!="",]
circ_data <- circ_data[circ_data$ctDNA.postTx=="POSITIVE",]
median_ctDNA <- median(circ_data$ctDNA.postTx.MTM, na.rm = TRUE)
range_ctDNA <- range(circ_data$ctDNA.postTx.MTM, na.rm = TRUE)
cat("Median MTM/mL post-treatment:", median_ctDNA, "\n")
Median MTM/mL post-treatment: 7.9
cat("Range MTM/mL post-treatment:", range_ctDNA[1], "-", range_ctDNA[2], "\n")
Range MTM/mL post-treatment: 0.1 - 737.8
#Median time from end treatment to progression for ctDNA negative pts post-treatment
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.postTx!="",]
circ_data <- circ_data[circ_data$ctDNA.postTx=="NEGATIVE",]
circ_data <- circ_data[circ_data$PFS.Event=="TRUE",]
median_PFS <- median(circ_data$PFS.months, na.rm = TRUE)
range_PFS <- range(circ_data$PFS.months, na.rm = TRUE)
cat("Median PFS:", median_PFS, "\n")
Median PFS: 4.977494
cat("Range PFS:", range_PFS[1], "-", range_PFS[2], "\n")
Range PFS: 3.055492 - 19.67999
#PFS by ctDNA status post/during-ICI
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.postTx!="",]
circ_datadf <- as.data.frame(circ_data)
survfit(Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event)~ctDNA.postTx, data = circ_data)
Call: survfit(formula = Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event) ~
ctDNA.postTx, data = circ_data)
n events median 0.95LCL 0.95UCL
ctDNA.postTx=NEGATIVE 9 4 NA 6.70 NA
ctDNA.postTx=POSITIVE 20 19 3.81 2.37 6.64
event_summary <- circ_data %>%
group_by(ctDNA.postTx) %>%
summarise(
Total = n(),
Events = sum(PFS.Event),
Fraction = Events / n(),
Percentage = (Events / n()) * 100
)
print(event_summary)
surv_object <-Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.postTx, data = circ_data,conf.int=0.95,conf.type="log-log")
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","red"), title="PFS - ctDNA Post/Under treatment", ylab= "Progression-Free Survival", xlab="Months from Start of Immunotherapy", legend.labs=c("ctDNA Negative", "ctDNA Positive"), legend.title="")
summary(KM_curve, times= c(0, 12, 24, 36))
Call: survfit(formula = surv_object ~ ctDNA.postTx, data = circ_data,
conf.int = 0.95, conf.type = "log-log")
ctDNA.postTx=NEGATIVE
time n.risk n.event survival std.err lower 95% CI upper 95% CI
0 9 0 1.000 0.000 1.000 1.000
12 5 3 0.667 0.157 0.282 0.878
24 3 1 0.533 0.173 0.177 0.796
36 2 0 0.533 0.173 0.177 0.796
ctDNA.postTx=POSITIVE
time n.risk n.event survival std.err lower 95% CI upper 95% CI
0 20 0 1.00 0.0000 1.00000 1.000
12 2 18 0.10 0.0671 0.01698 0.272
24 1 1 0.05 0.0487 0.00345 0.205
36 1 0 0.05 0.0487 0.00345 0.205
circ_data$ctDNA.postTx <- factor(circ_data$ctDNA.postTx, levels=c("NEGATIVE","POSITIVE"))
cox_fit <- coxph(surv_object ~ ctDNA.postTx, data=circ_data)
ggforest(cox_fit,data = circ_data)
summary(cox_fit)
Call:
coxph(formula = surv_object ~ ctDNA.postTx, data = circ_data)
n= 29, number of events= 23
coef exp(coef) se(coef) z Pr(>|z|)
ctDNA.postTxPOSITIVE 1.5616 4.7666 0.5703 2.738 0.00618 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95
ctDNA.postTxPOSITIVE 4.767 0.2098 1.559 14.58
Concordance= 0.649 (se = 0.049 )
Likelihood ratio test= 9.7 on 1 df, p=0.002
Wald test = 7.5 on 1 df, p=0.006
Score (logrank) test = 8.78 on 1 df, p=0.003
cox_fit_summary <- summary(cox_fit)
#Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)
[1] "HR = 4.77 (1.56-14.58); p = 0.006"
circ_data$ctDNA.postTx <- factor(circ_data$ctDNA.postTx, levels = c("NEGATIVE", "POSITIVE"), labels = c("Negative", "Positive"))
circ_data$PFS.Event <- factor(circ_data$PFS.Event, levels = c("FALSE", "TRUE"), labels = c("No Progression", "Progression"))
contingency_table <- table(circ_data$ctDNA.postTx, circ_data$PFS.Event)
chi_square_test <- chisq.test(contingency_table)
Warning: Chi-squared approximation may be incorrect
print(chi_square_test)
Pearson's Chi-squared test with Yates' continuity correction
data: contingency_table
X-squared = 6.8323, df = 1, p-value = 0.008952
fisher_exact_test <- fisher.test(contingency_table)
print(fisher_exact_test)
Fisher's Exact Test for Count Data
data: contingency_table
p-value = 0.005482
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
1.68957 1169.09723
sample estimates:
odds ratio
20.33413
print(contingency_table)
No Progression Progression
Negative 5 4
Positive 1 19
table_df <- as.data.frame(contingency_table)
table_df$Total <- ave(table_df$Freq, table_df$Var1, FUN = sum)
table_df$Percentage <- table_df$Freq / table_df$Total
table_df$MiddlePercentage <- table_df$Percentage / 2
ggplot(table_df, aes(x = Var1, y = Percentage, fill = Var2)) +
geom_bar(stat = "identity") +
geom_text(aes(y = MiddlePercentage, label = Freq), position = "stack", color = "black", vjust = 1.5, size = 7) +
theme_minimal() +
labs(title = "ctDNA status post/under-treatment",
x = "ctDNA",
y = "Patients (%)",
fill = "Progression",
caption = paste("Fisher's exact test p-value: ", format.pval(fisher_exact_test$p.value))) +
scale_y_continuous(labels = scales::percent_format()) +
scale_fill_manual(values = c("No Progression" = "blue", "Progression" = "red")) + # define custom colors
theme(axis.text.x = element_text(angle = 0, hjust = 1.5, size = 14), # increase x-axis text size
axis.text.y = element_text(size = 14, color = "black"), # increase y-axis text size
axis.title.x = element_text(size = 14, color = "black"), # increase x-axis label size
axis.title.y = element_text(size = 14, color = "black"), # increase y-axis label size
legend.text = element_text(size = 12, color = "black")) # increase Progression label size
#OS by ctDNA status post/during-ICI
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.postTx!="",]
circ_datadf <- as.data.frame(circ_data)
survfit(Surv(time = circ_data$OS.months, event = circ_data$OS.Event)~ctDNA.postTx, data = circ_data)
Call: survfit(formula = Surv(time = circ_data$OS.months, event = circ_data$OS.Event) ~
ctDNA.postTx, data = circ_data)
n events median 0.95LCL 0.95UCL
ctDNA.postTx=NEGATIVE 9 4 NA 16.85 NA
ctDNA.postTx=POSITIVE 20 17 10.4 7.75 18
event_summary <- circ_data %>%
group_by(ctDNA.postTx) %>%
summarise(
Total = n(),
Events = sum(OS.Event),
Fraction = Events / n(),
Percentage = (Events / n()) * 100
)
print(event_summary)
surv_object <-Surv(time = circ_data$OS.months, event = circ_data$OS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.postTx, data = circ_data,conf.int=0.95,conf.type="log-log")
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","red"), title="OS - ctDNA Post/Under treatment", ylab= "Overall Survival", xlab="Months from Start of Immunotherapy", legend.labs=c("ctDNA Negative", "ctDNA Positive"), legend.title="")
summary(KM_curve, times= c(12, 24, 36))
Call: survfit(formula = surv_object ~ ctDNA.postTx, data = circ_data,
conf.int = 0.95, conf.type = "log-log")
ctDNA.postTx=NEGATIVE
time n.risk n.event survival std.err lower 95% CI upper 95% CI
12 7 1 0.889 0.105 0.433 0.984
24 3 3 0.508 0.177 0.157 0.781
36 2 0 0.508 0.177 0.157 0.781
ctDNA.postTx=POSITIVE
time n.risk n.event survival std.err lower 95% CI upper 95% CI
12 9 11 0.45 0.1112 0.2311 0.647
24 3 5 0.18 0.0900 0.0480 0.380
36 2 1 0.12 0.0775 0.0213 0.311
circ_data$ctDNA.postTx <- factor(circ_data$ctDNA.postTx, levels=c("NEGATIVE","POSITIVE"))
cox_fit <- coxph(surv_object ~ ctDNA.postTx, data=circ_data)
ggforest(cox_fit,data = circ_data)
summary(cox_fit)
Call:
coxph(formula = surv_object ~ ctDNA.postTx, data = circ_data)
n= 29, number of events= 21
coef exp(coef) se(coef) z Pr(>|z|)
ctDNA.postTxPOSITIVE 1.2213 3.3918 0.5621 2.173 0.0298 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95
ctDNA.postTxPOSITIVE 3.392 0.2948 1.127 10.21
Concordance= 0.638 (se = 0.048 )
Likelihood ratio test= 5.81 on 1 df, p=0.02
Wald test = 4.72 on 1 df, p=0.03
Score (logrank) test = 5.29 on 1 df, p=0.02
cox_fit_summary <- summary(cox_fit)
#Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)
[1] "HR = 3.39 (1.13-10.21); p = 0.03"
circ_data$ctDNA.postTx <- factor(circ_data$ctDNA.postTx, levels = c("NEGATIVE", "POSITIVE"), labels = c("Negative", "Positive"))
circ_data$OS.Event <- factor(circ_data$OS.Event, levels = c("FALSE", "TRUE"), labels = c("Alive", "Deceased"))
contingency_table <- table(circ_data$ctDNA.postTx, circ_data$OS.Event)
chi_square_test <- chisq.test(contingency_table)
Warning: Chi-squared approximation may be incorrect
print(chi_square_test)
Pearson's Chi-squared test with Yates' continuity correction
data: contingency_table
X-squared = 3.2819, df = 1, p-value = 0.07005
fisher_exact_test <- fisher.test(contingency_table)
print(fisher_exact_test)
Fisher's Exact Test for Count Data
data: contingency_table
p-value = 0.0667
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
0.8624077 62.5051555
sample estimates:
odds ratio
6.505452
print(contingency_table)
Alive Deceased
Negative 5 4
Positive 3 17
table_df <- as.data.frame(contingency_table)
table_df$Total <- ave(table_df$Freq, table_df$Var1, FUN = sum)
table_df$Percentage <- table_df$Freq / table_df$Total
table_df$MiddlePercentage <- table_df$Percentage / 2
ggplot(table_df, aes(x = Var1, y = Percentage, fill = Var2)) +
geom_bar(stat = "identity") +
geom_text(aes(y = MiddlePercentage, label = Freq), position = "stack", color = "black", vjust = 1.5, size = 7) +
theme_minimal() +
labs(title = "ctDNA status post/under-treatment",
x = "ctDNA",
y = "Patients (%)",
fill = "Living Status",
caption = paste("Fisher's exact test p-value: ", format.pval(fisher_exact_test$p.value))) +
scale_y_continuous(labels = scales::percent_format()) +
scale_fill_manual(values = c("Alive" = "blue", "Deceased" = "red")) + # define custom colors
theme(axis.text.x = element_text(angle = 0, hjust = 1.5, size = 14), # increase x-axis text size
axis.text.y = element_text(size = 14, color = "black"), # increase y-axis text size
axis.title.x = element_text(size = 14, color = "black"), # increase x-axis label size
axis.title.y = element_text(size = 14, color = "black"), # increase y-axis label size
legend.text = element_text(size = 12, color = "black")) # increase Progression label size
#Association of ctDNA status post/during-ICI with BOR
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_datadf <- as.data.frame(circ_data)
circ_data$ctDNA.postTx <- factor(circ_data$ctDNA.postTx, levels = c("NEGATIVE", "POSITIVE"), labels = c("Negative", "Positive"))
circ_data$RESIST <- factor(circ_data$RESIST, levels = c("CR", "PR", "SD", "PD"))
contingency_table <- table(circ_data$ctDNA.postTx, circ_data$RESIST)
chi_square_test <- chisq.test(contingency_table)
Warning: Chi-squared approximation may be incorrect
print(chi_square_test)
Pearson's Chi-squared test
data: contingency_table
X-squared = 11.869, df = 3, p-value = 0.007847
fisher_exact_test <- fisher.test(contingency_table)
print(fisher_exact_test)
Fisher's Exact Test for Count Data
data: contingency_table
p-value = 0.006541
alternative hypothesis: two.sided
print(contingency_table)
CR PR SD PD
Negative 4 3 0 2
Positive 1 2 7 10
table_df <- as.data.frame(contingency_table)
table_df$Total <- ave(table_df$Freq, table_df$Var1, FUN = sum)
table_df$Percentage <- table_df$Freq / table_df$Total
table_df$MiddlePercentage <- table_df$Percentage / 2
ggplot(table_df, aes(x = Var1, y = Percentage, fill = Var2)) +
geom_bar(stat = "identity") +
geom_text(aes(y = MiddlePercentage, label = Freq), position = "stack", color = "black", vjust = 1.5, size = 7) +
theme_minimal() +
labs(title = "ctDNA status post/under-treatment",
x = "ctDNA",
y = "Patients (%)",
fill = "BOR",
caption = paste("Fisher's exact test p-value: ", format.pval(fisher_exact_test$p.value))) +
scale_y_continuous(labels = scales::percent_format()) +
scale_fill_manual(values = c("CR" = "blue", "PR" = "lightblue", "SD" = "orange", "PD" = "red")) + # define custom colors
theme(axis.text.x = element_text(angle = 0, hjust = 1.5, size = 14), # increase x-axis text size
axis.text.y = element_text(size = 14, color = "black"), # increase y-axis text size
axis.title.x = element_text(size = 14, color = "black"), # increase x-axis label size
axis.title.y = element_text(size = 14, color = "black"), # increase y-axis label size
legend.text = element_text(size = 12, color = "black")) # increase Progression label size
#PFS by ctDNA kinetics post/during-ICI
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ΔctDNA!="",]
circ_datadf <- as.data.frame(circ_data)
survfit(Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event)~ΔctDNA, data = circ_data)
Call: survfit(formula = Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event) ~
ΔctDNA, data = circ_data)
n events median 0.95LCL 0.95UCL
ΔctDNA=NEGATIVE 13 8 18.04 6.7 NA
ΔctDNA=POSITIVE 12 12 2.41 2.3 NA
event_summary <- circ_data %>%
group_by(ΔctDNA) %>%
summarise(
Total = n(),
Events = sum(PFS.Event),
Fraction = Events / n(),
Percentage = (Events / n()) * 100
)
print(event_summary)
surv_object <-Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event)
KM_curve <- survfit(surv_object ~ ΔctDNA, data = circ_data,conf.int=0.95,conf.type="log-log")
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","red"), title="PFS - ctDNA Kinetics Post/Under treatment", ylab= "Progression-Free Survival", xlab="Months from Start of Immunotherapy", legend.labs=c("Decreased ΔctDNA", "Increased ΔctDNA"), legend.title="")
summary(KM_curve, times= c(0, 12, 24, 36))
Call: survfit(formula = surv_object ~ ΔctDNA, data = circ_data, conf.int = 0.95,
conf.type = "log-log")
ΔctDNA=NEGATIVE
time n.risk n.event survival std.err lower 95% CI upper 95% CI
0 13 0 1.000 0.000 1.000 1.000
12 6 6 0.538 0.138 0.248 0.760
24 3 2 0.359 0.139 0.117 0.613
36 2 0 0.359 0.139 0.117 0.613
ΔctDNA=POSITIVE
time n.risk n.event survival std.err lower 95% CI upper 95% CI
0 12 0 1 0 1 1
circ_data$ΔctDNA <- factor(circ_data$ΔctDNA, levels=c("NEGATIVE","POSITIVE"))
cox_fit <- coxph(surv_object ~ ΔctDNA, data=circ_data)
ggforest(cox_fit,data = circ_data)
summary(cox_fit)
Call:
coxph(formula = surv_object ~ ΔctDNA, data = circ_data)
n= 25, number of events= 20
coef exp(coef) se(coef) z Pr(>|z|)
ΔctDNAPOSITIVE 2.716 15.112 0.794 3.42 0.000626 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95
ΔctDNAPOSITIVE 15.11 0.06617 3.188 71.64
Concordance= 0.716 (se = 0.044 )
Likelihood ratio test= 17.9 on 1 df, p=2e-05
Wald test = 11.7 on 1 df, p=6e-04
Score (logrank) test = 18.33 on 1 df, p=2e-05
cox_fit_summary <- summary(cox_fit)
#Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)
[1] "HR = 15.11 (3.19-71.64); p = 0.001"
circ_data$ΔctDNA <- factor(circ_data$ΔctDNA, levels = c("NEGATIVE", "POSITIVE"), labels = c("Decreased", "Increased"))
circ_data$PFS.Event <- factor(circ_data$PFS.Event, levels = c("FALSE", "TRUE"), labels = c("No Progression", "Progression"))
contingency_table <- table(circ_data$ΔctDNA, circ_data$PFS.Event)
chi_square_test <- chisq.test(contingency_table)
Warning: Chi-squared approximation may be incorrect
print(chi_square_test)
Pearson's Chi-squared test with Yates' continuity correction
data: contingency_table
X-squared = 3.6158, df = 1, p-value = 0.05723
fisher_exact_test <- fisher.test(contingency_table)
print(fisher_exact_test)
Fisher's Exact Test for Count Data
data: contingency_table
p-value = 0.03913
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
1.013281 Inf
sample estimates:
odds ratio
Inf
print(contingency_table)
No Progression Progression
Decreased 5 8
Increased 0 12
table_df <- as.data.frame(contingency_table)
table_df$Total <- ave(table_df$Freq, table_df$Var1, FUN = sum)
table_df$Percentage <- table_df$Freq / table_df$Total
table_df$MiddlePercentage <- table_df$Percentage / 2
ggplot(table_df, aes(x = Var1, y = Percentage, fill = Var2)) +
geom_bar(stat = "identity") +
geom_text(aes(y = MiddlePercentage, label = Freq), position = "stack", color = "black", vjust = 1.5, size = 7) +
theme_minimal() +
labs(title = "ctDNA kinetics post/under-treatment",
x = "ctDNA",
y = "Patients (%)",
fill = "Progression",
caption = paste("Fisher's exact test p-value: ", format.pval(fisher_exact_test$p.value))) +
scale_y_continuous(labels = scales::percent_format()) +
scale_fill_manual(values = c("No Progression" = "blue", "Progression" = "red")) + # define custom colors
theme(axis.text.x = element_text(angle = 0, hjust = 1.5, size = 14), # increase x-axis text size
axis.text.y = element_text(size = 14, color = "black"), # increase y-axis text size
axis.title.x = element_text(size = 14, color = "black"), # increase x-axis label size
axis.title.y = element_text(size = 14, color = "black"), # increase y-axis label size
legend.text = element_text(size = 12, color = "black")) # increase Progression label size
#OS by ctDNA kinetics post/during-ICI
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ΔctDNA!="",]
circ_datadf <- as.data.frame(circ_data)
survfit(Surv(time = circ_data$OS.months, event = circ_data$OS.Event)~ΔctDNA, data = circ_data)
Call: survfit(formula = Surv(time = circ_data$OS.months, event = circ_data$OS.Event) ~
ΔctDNA, data = circ_data)
n events median 0.95LCL 0.95UCL
ΔctDNA=NEGATIVE 13 8 18.0 9.82 NA
ΔctDNA=POSITIVE 12 10 10.8 4.44 NA
event_summary <- circ_data %>%
group_by(ΔctDNA) %>%
summarise(
Total = n(),
Events = sum(OS.Event),
Fraction = Events / n(),
Percentage = (Events / n()) * 100
)
print(event_summary)
surv_object <-Surv(time = circ_data$OS.months, event = circ_data$OS.Event)
KM_curve <- survfit(surv_object ~ ΔctDNA, data = circ_data,conf.int=0.95,conf.type="log-log")
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","red"), title="OS - ctDNA Kinetics Post/Under treatment", ylab= "Overall Survival", xlab="Months from Start of Immunotherapy", legend.labs=c("Decreased ΔctDNA", "Increased ΔctDNA"), legend.title="")
summary(KM_curve, times= c(0, 12, 24, 36))
Call: survfit(formula = surv_object ~ ΔctDNA, data = circ_data, conf.int = 0.95,
conf.type = "log-log")
ΔctDNA=NEGATIVE
time n.risk n.event survival std.err lower 95% CI upper 95% CI
0 13 0 1.000 0.000 1.000 1.000
12 7 5 0.615 0.135 0.308 0.818
24 3 3 0.352 0.139 0.112 0.607
36 2 0 0.352 0.139 0.112 0.607
ΔctDNA=POSITIVE
time n.risk n.event survival std.err lower 95% CI upper 95% CI
0 12 0 1.000 0.000 1.00000 1.000
12 6 6 0.500 0.144 0.20848 0.736
24 2 3 0.222 0.128 0.04111 0.492
36 1 1 0.111 0.101 0.00701 0.378
circ_data$ΔctDNA <- factor(circ_data$ΔctDNA, levels=c("NEGATIVE","POSITIVE"))
cox_fit <- coxph(surv_object ~ ΔctDNA, data=circ_data)
ggforest(cox_fit,data = circ_data)
summary(cox_fit)
Call:
coxph(formula = surv_object ~ ΔctDNA, data = circ_data)
n= 25, number of events= 18
coef exp(coef) se(coef) z Pr(>|z|)
ΔctDNAPOSITIVE 0.6947 2.0031 0.4806 1.445 0.148
exp(coef) exp(-coef) lower .95 upper .95
ΔctDNAPOSITIVE 2.003 0.4992 0.7809 5.138
Concordance= 0.598 (se = 0.061 )
Likelihood ratio test= 2.1 on 1 df, p=0.1
Wald test = 2.09 on 1 df, p=0.1
Score (logrank) test = 2.17 on 1 df, p=0.1
cox_fit_summary <- summary(cox_fit)
#Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)
[1] "HR = 2 (0.78-5.14); p = 0.148"
circ_data$ΔctDNA <- factor(circ_data$ΔctDNA, levels = c("NEGATIVE", "POSITIVE"), labels = c("Decreased", "Increased"))
circ_data$OS.Event <- factor(circ_data$OS.Event, levels = c("FALSE", "TRUE"), labels = c("Alive", "Deceased"))
contingency_table <- table(circ_data$ΔctDNA, circ_data$OS.Event)
chi_square_test <- chisq.test(contingency_table)
Warning: Chi-squared approximation may be incorrect
print(chi_square_test)
Pearson's Chi-squared test with Yates' continuity correction
data: contingency_table
X-squared = 0.58792, df = 1, p-value = 0.4432
fisher_exact_test <- fisher.test(contingency_table)
print(fisher_exact_test)
Fisher's Exact Test for Count Data
data: contingency_table
p-value = 0.3783
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
0.363132 39.390335
sample estimates:
odds ratio
2.985191
print(contingency_table)
Alive Deceased
Decreased 5 8
Increased 2 10
table_df <- as.data.frame(contingency_table)
table_df$Total <- ave(table_df$Freq, table_df$Var1, FUN = sum)
table_df$Percentage <- table_df$Freq / table_df$Total
table_df$MiddlePercentage <- table_df$Percentage / 2
ggplot(table_df, aes(x = Var1, y = Percentage, fill = Var2)) +
geom_bar(stat = "identity") +
geom_text(aes(y = MiddlePercentage, label = Freq), position = "stack", color = "black", vjust = 1.5, size = 7) +
theme_minimal() +
labs(title = "ctDNA kinetics post/under-treatment",
x = "ctDNA",
y = "Patients (%)",
fill = "Vital Status",
caption = paste("Fisher's exact test p-value: ", format.pval(fisher_exact_test$p.value))) +
scale_y_continuous(labels = scales::percent_format()) +
scale_fill_manual(values = c("Alive" = "blue", "Deceased" = "red")) + # define custom colors
theme(axis.text.x = element_text(angle = 0, hjust = 1.5, size = 14), # increase x-axis text size
axis.text.y = element_text(size = 14, color = "black"), # increase y-axis text size
axis.title.x = element_text(size = 14, color = "black"), # increase x-axis label size
axis.title.y = element_text(size = 14, color = "black"), # increase y-axis label size
legend.text = element_text(size = 12, color = "black")) # increase Progression label size
#Association of ctDNA kinetics post/during-ICI with BOR
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ΔctDNA!="",]
circ_datadf <- as.data.frame(circ_data)
circ_data$ΔctDNA <- factor(circ_data$ΔctDNA, levels = c("NEGATIVE", "POSITIVE"), labels = c("Decreased", "Increased"))
circ_data$RESIST <- factor(circ_data$RESIST, levels = c("CR", "PR", "SD", "PD"))
contingency_table <- table(circ_data$ΔctDNA, circ_data$RESIST)
chi_square_test <- chisq.test(contingency_table)
Warning: Chi-squared approximation may be incorrect
print(chi_square_test)
Pearson's Chi-squared test
data: contingency_table
X-squared = 15.385, df = 3, p-value = 0.001516
fisher_exact_test <- fisher.test(contingency_table)
print(fisher_exact_test)
Fisher's Exact Test for Count Data
data: contingency_table
p-value = 0.0003154
alternative hypothesis: two.sided
print(contingency_table)
CR PR SD PD
Decreased 4 5 3 1
Increased 0 0 3 9
table_df <- as.data.frame(contingency_table)
table_df$Total <- ave(table_df$Freq, table_df$Var1, FUN = sum)
table_df$Percentage <- table_df$Freq / table_df$Total
table_df$MiddlePercentage <- table_df$Percentage / 2
ggplot(table_df, aes(x = Var1, y = Percentage, fill = Var2)) +
geom_bar(stat = "identity") +
geom_text(aes(y = MiddlePercentage, label = Freq), position = "stack", color = "black", vjust = 1.5, size = 7) +
theme_minimal() +
labs(title = "ctDNA kinetics post/under-treatment",
x = "ctDNA",
y = "Patients (%)",
fill = "BOR",
caption = paste("Fisher's exact test p-value: ", format.pval(fisher_exact_test$p.value))) +
scale_y_continuous(labels = scales::percent_format()) +
scale_fill_manual(values = c("CR" = "blue", "PR" = "lightblue", "SD" = "orange", "PD" = "red")) + # define custom colors
theme(axis.text.x = element_text(angle = 0, hjust = 1.5, size = 14), # increase x-axis text size
axis.text.y = element_text(size = 14, color = "black"), # increase y-axis text size
axis.title.x = element_text(size = 14, color = "black"), # increase x-axis label size
axis.title.y = element_text(size = 14, color = "black"), # increase y-axis label size
legend.text = element_text(size = 12, color = "black")) # increase Progression label size
#PFS by ctDNA clearance post/during-ICI
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Base!="",]
circ_data <- circ_data[circ_data$ctDNA.Base=="POSITIVE",]
circ_datadf <- as.data.frame(circ_data)
circ_data$ctDNA.Dynamics <- NA #first we create the variable for the ctDNA & NAC combination, and we assign values
circ_data <- circ_data %>%
mutate(ctDNA.Dynamics = case_when(
ctDNA.Base == "POSITIVE" & ctDNA.postTx == "NEGATIVE" ~ 1,
ctDNA.Base == "POSITIVE" & ctDNA.postTx == "POSITIVE" ~ 2
))
survfit(Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event)~ctDNA.Dynamics, data = circ_data)
Call: survfit(formula = Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event) ~
ctDNA.Dynamics, data = circ_data)
n events median 0.95LCL 0.95UCL
ctDNA.Dynamics=1 6 2 NA 19.68 NA
ctDNA.Dynamics=2 19 18 3.68 2.37 7.95
event_summary <- circ_data %>%
group_by(ctDNA.Dynamics) %>%
summarise(
Total = n(),
Events = sum(PFS.Event),
Fraction = Events / n(),
Percentage = (Events / n()) * 100
)
print(event_summary)
surv_object <-Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.Dynamics, data = circ_data,conf.int=0.95,conf.type="log-log")
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","red"), title="PFS - ctDNA clearance Post/Under treatment", ylab= "Progression-Free Survival", xlab="Months from Start of Immunotherapy", legend.labs=c("Clearance", "No Clearance"), legend.title="")
summary(KM_curve, times= c(0, 12, 24, 36))
Call: survfit(formula = surv_object ~ ctDNA.Dynamics, data = circ_data,
conf.int = 0.95, conf.type = "log-log")
ctDNA.Dynamics=1
time n.risk n.event survival std.err lower 95% CI upper 95% CI
0 6 0 1.000 0.000 1.000 1.000
12 4 1 0.833 0.152 0.273 0.975
24 2 1 0.625 0.213 0.142 0.893
36 1 0 0.625 0.213 0.142 0.893
ctDNA.Dynamics=2
time n.risk n.event survival std.err lower 95% CI upper 95% CI
0 19 0 1.0000 0.0000 1.00000 1.000
12 2 17 0.1053 0.0704 0.01777 0.284
24 1 1 0.0526 0.0512 0.00359 0.214
36 1 0 0.0526 0.0512 0.00359 0.214
circ_data$ctDNA.Dynamics <- factor(circ_data$ctDNA.Dynamics, levels=c("1","2"), labels = c("Clearance", "No Clearance"))
cox_fit <- coxph(surv_object ~ ctDNA.Dynamics, data=circ_data)
ggforest(cox_fit,data = circ_data)
summary(cox_fit)
Call:
coxph(formula = surv_object ~ ctDNA.Dynamics, data = circ_data)
n= 25, number of events= 20
coef exp(coef) se(coef) z Pr(>|z|)
ctDNA.DynamicsNo Clearance 2.085 8.048 0.767 2.719 0.00655 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95
ctDNA.DynamicsNo Clearance 8.048 0.1243 1.79 36.19
Concordance= 0.672 (se = 0.049 )
Likelihood ratio test= 11.61 on 1 df, p=7e-04
Wald test = 7.39 on 1 df, p=0.007
Score (logrank) test = 9.78 on 1 df, p=0.002
cox_fit_summary <- summary(cox_fit)
# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)
[1] "HR = 8.05 (1.79-36.19); p = 0.007"
circ_data$PFS.Event <- factor(circ_data$PFS.Event, levels = c("FALSE", "TRUE"), labels = c("No Progression", "Progression"))
contingency_table <- table(circ_data$ctDNA.Dynamics, circ_data$PFS.Event)
chi_square_test <- chisq.test(contingency_table)
Warning: Chi-squared approximation may be incorrect
print(chi_square_test)
Pearson's Chi-squared test with Yates' continuity correction
data: contingency_table
X-squared = 7.2505, df = 1, p-value = 0.007088
fisher_exact_test <- fisher.test(contingency_table)
print(fisher_exact_test)
Fisher's Exact Test for Count Data
data: contingency_table
p-value = 0.005477
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
1.796359 1833.444857
sample estimates:
odds ratio
27.59073
print(contingency_table)
No Progression Progression
Clearance 4 2
No Clearance 1 18
table_df <- as.data.frame(contingency_table)
table_df$Total <- ave(table_df$Freq, table_df$Var1, FUN = sum)
table_df$Percentage <- table_df$Freq / table_df$Total
table_df$MiddlePercentage <- table_df$Percentage / 2
ggplot(table_df, aes(x = Var1, y = Percentage, fill = Var2)) +
geom_bar(stat = "identity") +
geom_text(aes(y = MiddlePercentage, label = Freq), position = "stack", color = "black", vjust = 1.5, size = 7) +
theme_minimal() +
labs(title = "ctDNA clearance post/under-treatment",
x = "ctDNA",
y = "Patients (%)",
fill = "Progression",
caption = paste("Fisher's exact test p-value: ", format.pval(fisher_exact_test$p.value))) +
scale_y_continuous(labels = scales::percent_format()) +
scale_fill_manual(values = c("No Progression" = "blue", "Progression" = "red")) + # define custom colors
theme(axis.text.x = element_text(angle = 0, hjust = 1.5, size = 14), # increase x-axis text size
axis.text.y = element_text(size = 14, color = "black"), # increase y-axis text size
axis.title.x = element_text(size = 14, color = "black"), # increase x-axis label size
axis.title.y = element_text(size = 14, color = "black"), # increase y-axis label size
legend.text = element_text(size = 12, color = "black")) # increase Progression label size
#OS by ctDNA clearance post/during-ICI
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Base!="",]
circ_data <- circ_data[circ_data$ctDNA.Base=="POSITIVE",]
circ_datadf <- as.data.frame(circ_data)
circ_data$ctDNA.Dynamics <- NA #first we create the variable for the ctDNA & NAC combination, and we assign values
circ_data <- circ_data %>%
mutate(ctDNA.Dynamics = case_when(
ctDNA.Base == "POSITIVE" & ctDNA.postTx == "NEGATIVE" ~ 1,
ctDNA.Base == "POSITIVE" & ctDNA.postTx == "POSITIVE" ~ 2
))
survfit(Surv(time = circ_data$OS.months, event = circ_data$OS.Event)~ctDNA.Dynamics, data = circ_data)
Call: survfit(formula = Surv(time = circ_data$OS.months, event = circ_data$OS.Event) ~
ctDNA.Dynamics, data = circ_data)
n events median 0.95LCL 0.95UCL
ctDNA.Dynamics=1 6 2 NA 20.50 NA
ctDNA.Dynamics=2 19 16 9.82 7.75 26.7
event_summary <- circ_data %>%
group_by(ctDNA.Dynamics) %>%
summarise(
Total = n(),
Events = sum(OS.Event),
Fraction = Events / n(),
Percentage = (Events / n()) * 100
)
print(event_summary)
surv_object <-Surv(time = circ_data$OS.months, event = circ_data$OS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.Dynamics, data = circ_data,conf.int=0.95,conf.type="log-log")
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","red"), title="OS - ctDNA clearance Post/Under treatment", ylab= "Overall Survival", xlab="Months from Start of Immunotherapy", legend.labs=c("Clearance", "No Clearance"), legend.title="")
summary(KM_curve, times= c(0, 12, 24, 36))
Call: survfit(formula = surv_object ~ ctDNA.Dynamics, data = circ_data,
conf.int = 0.95, conf.type = "log-log")
ctDNA.Dynamics=1
time n.risk n.event survival std.err lower 95% CI upper 95% CI
0 6 0 1.0 0.000 1.000 1.000
12 5 0 1.0 0.000 NA NA
24 2 2 0.6 0.219 0.126 0.882
36 1 0 0.6 0.219 0.126 0.882
ctDNA.Dynamics=2
time n.risk n.event survival std.err lower 95% CI upper 95% CI
0 19 0 1.000 0.0000 1.0000 1.000
12 8 11 0.421 0.1133 0.2037 0.625
24 3 4 0.189 0.0942 0.0503 0.396
36 2 1 0.126 0.0813 0.0222 0.325
circ_data$ctDNA.Dynamics <- factor(circ_data$ctDNA.Dynamics, levels=c("1","2"), labels = c("Clearance", "No Clearance"))
cox_fit <- coxph(surv_object ~ ctDNA.Dynamics, data=circ_data)
ggforest(cox_fit,data = circ_data)
summary(cox_fit)
Call:
coxph(formula = surv_object ~ ctDNA.Dynamics, data = circ_data)
n= 25, number of events= 18
coef exp(coef) se(coef) z Pr(>|z|)
ctDNA.DynamicsNo Clearance 1.6133 5.0196 0.7578 2.129 0.0333 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95
ctDNA.DynamicsNo Clearance 5.02 0.1992 1.137 22.17
Concordance= 0.65 (se = 0.047 )
Likelihood ratio test= 6.57 on 1 df, p=0.01
Wald test = 4.53 on 1 df, p=0.03
Score (logrank) test = 5.52 on 1 df, p=0.02
cox_fit_summary <- summary(cox_fit)
# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)
[1] "HR = 5.02 (1.14-22.17); p = 0.033"
circ_data$OS.Event <- factor(circ_data$OS.Event, levels = c("FALSE", "TRUE"), labels = c("Alive", "Deceased"))
contingency_table <- table(circ_data$ctDNA.Dynamics, circ_data$OS.Event)
chi_square_test <- chisq.test(contingency_table)
Warning: Chi-squared approximation may be incorrect
print(chi_square_test)
Pearson's Chi-squared test with Yates' continuity correction
data: contingency_table
X-squared = 3.6032, df = 1, p-value = 0.05767
fisher_exact_test <- fisher.test(contingency_table)
print(fisher_exact_test)
Fisher's Exact Test for Count Data
data: contingency_table
p-value = 0.03241
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
0.8988842 150.9697807
sample estimates:
odds ratio
9.34674
print(contingency_table)
Alive Deceased
Clearance 4 2
No Clearance 3 16
table_df <- as.data.frame(contingency_table)
table_df$Total <- ave(table_df$Freq, table_df$Var1, FUN = sum)
table_df$Percentage <- table_df$Freq / table_df$Total
table_df$MiddlePercentage <- table_df$Percentage / 2
ggplot(table_df, aes(x = Var1, y = Percentage, fill = Var2)) +
geom_bar(stat = "identity") +
geom_text(aes(y = MiddlePercentage, label = Freq), position = "stack", color = "black", vjust = 1.5, size = 7) +
theme_minimal() +
labs(title = "ctDNA clearance post/under-treatment",
x = "ctDNA",
y = "Patients (%)",
fill = "Vital Status",
caption = paste("Fisher's exact test p-value: ", format.pval(fisher_exact_test$p.value))) +
scale_y_continuous(labels = scales::percent_format()) +
scale_fill_manual(values = c("Alive" = "blue", "Deceased" = "red")) + # define custom colors
theme(axis.text.x = element_text(angle = 0, hjust = 1.5, size = 14), # increase x-axis text size
axis.text.y = element_text(size = 14, color = "black"), # increase y-axis text size
axis.title.x = element_text(size = 14, color = "black"), # increase x-axis label size
axis.title.y = element_text(size = 14, color = "black"), # increase y-axis label size
legend.text = element_text(size = 12, color = "black")) # increase Progression label size
#Association of ctDNA clearance post/during-ICI with BOR
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Base!="",]
circ_data <- circ_data[circ_data$ctDNA.Base=="POSITIVE",]
circ_datadf <- as.data.frame(circ_data)
circ_data$ctDNA.Dynamics <- NA #first we create the variable for the ctDNA & NAC combination, and we assign values
circ_data <- circ_data %>%
mutate(ctDNA.Dynamics = case_when(
ctDNA.Base == "POSITIVE" & ctDNA.postTx == "NEGATIVE" ~ 1,
ctDNA.Base == "POSITIVE" & ctDNA.postTx == "POSITIVE" ~ 2
))
circ_data$ctDNA.Dynamics <- factor(circ_data$ctDNA.Dynamics, levels=c("1","2"), labels = c("Clearance", "No Clearance"))
circ_data$RESIST <- factor(circ_data$RESIST, levels = c("CR", "PR", "SD", "PD"))
contingency_table <- table(circ_data$ctDNA.Dynamics, circ_data$RESIST)
chi_square_test <- chisq.test(contingency_table)
Warning: Chi-squared approximation may be incorrect
print(chi_square_test)
Pearson's Chi-squared test
data: contingency_table
X-squared = 14.309, df = 3, p-value = 0.002513
fisher_exact_test <- fisher.test(contingency_table)
print(fisher_exact_test)
Fisher's Exact Test for Count Data
data: contingency_table
p-value = 0.001129
alternative hypothesis: two.sided
print(contingency_table)
CR PR SD PD
Clearance 3 3 0 0
No Clearance 1 2 6 10
table_df <- as.data.frame(contingency_table)
table_df$Total <- ave(table_df$Freq, table_df$Var1, FUN = sum)
table_df$Percentage <- table_df$Freq / table_df$Total
table_df$MiddlePercentage <- table_df$Percentage / 2
ggplot(table_df, aes(x = Var1, y = Percentage, fill = Var2)) +
geom_bar(stat = "identity") +
geom_text(aes(y = MiddlePercentage, label = Freq), position = "stack", color = "black", vjust = 1.5, size = 7) +
theme_minimal() +
labs(title = "ctDNA clearance post/under-treatment",
x = "ctDNA",
y = "Patients (%)",
fill = "BOR",
caption = paste("Fisher's exact test p-value: ", format.pval(fisher_exact_test$p.value))) +
scale_y_continuous(labels = scales::percent_format()) +
scale_fill_manual(values = c("CR" = "blue", "PR" = "lightblue", "SD" = "orange", "PD" = "red")) + # define custom colors
theme(axis.text.x = element_text(angle = 0, hjust = 1.5, size = 14), # increase x-axis text size
axis.text.y = element_text(size = 14, color = "black"), # increase y-axis text size
axis.title.x = element_text(size = 14, color = "black"), # increase x-axis label size
axis.title.y = element_text(size = 14, color = "black"), # increase y-axis label size
legend.text = element_text(size = 12, color = "black")) # increase Progression label size
#Multivariate cox regression for PFS
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Base!="",]
circ_data <- circ_data[circ_data$ctDNA.Base=="POSITIVE",]
circ_datadf <- as.data.frame(circ_data)
circ_data$ctDNA.Dynamics <- NA #first we create the variable for the ctDNA & NAC combination, and we assign values
circ_data <- circ_data %>%
mutate(ctDNA.Dynamics = case_when(
ctDNA.Base == "POSITIVE" & ctDNA.postTx == "NEGATIVE" ~ 1,
ctDNA.Base == "POSITIVE" & ctDNA.postTx == "POSITIVE" ~ 2
))
circ_data$ctDNA.Dynamics <- factor(circ_data$ctDNA.Dynamics, levels=c("1","2"), labels = c("Clearance", "No Clearance"))
circ_data$cStage <- factor(circ_data$cStage, levels = c("0-II", "III-IV"))
circ_data$CPS.Scorev2 <- factor(circ_data$CPS.Scorev2, levels = c("≥1", "<1"))
surv_object <- Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event)
cox_fit <- coxph(surv_object ~ ctDNA.Dynamics + Age + cStage + CPS.Scorev2, data=circ_data)
ggforest(cox_fit, data = circ_data, main = "Multivariate Regression Model for PFS", refLabel = "Reference Group")
test.ph <- cox.zph(cox_fit)
#Multivariate cox regression for PFS v2
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Base!="",]
circ_data <- circ_data[circ_data$ctDNA.Base=="POSITIVE",]
circ_datadf <- as.data.frame(circ_data)
circ_data$ctDNA.Dynamics <- NA #first we create the variable for the ctDNA & NAC combination, and we assign values
circ_data <- circ_data %>%
mutate(ctDNA.Dynamics = case_when(
ctDNA.Base == "POSITIVE" & ctDNA.postTx == "NEGATIVE" ~ 1,
ctDNA.Base == "POSITIVE" & ctDNA.postTx == "POSITIVE" ~ 2
))
circ_data$ctDNA.Dynamics <- factor(circ_data$ctDNA.Dynamics, levels=c("1","2"), labels = c("Clearance", "No Clearance"))
circ_data$Inclusion.status <- factor(circ_data$Inclusion.status, levels = c("Loco-regional", "Metastatic"))
circ_data$p16.status <- factor(circ_data$p16.status, levels = c("Negative", "Positive"))
circ_data$CPS.Scorev2 <- factor(circ_data$CPS.Scorev2, levels = c("≥1", "<1"))
surv_object <- Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event)
cox_fit <- coxph(surv_object ~ ctDNA.Dynamics + Age + Inclusion.status + p16.status + CPS.Scorev2, data=circ_data)
ggforest(cox_fit, data = circ_data, main = "Multivariate Regression Model for PFS", refLabel = "Reference Group")
test.ph <- cox.zph(cox_fit)
#Multivariate cox regression for PFS v3
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Base!="",]
circ_data <- circ_data[circ_data$ctDNA.Base=="POSITIVE",]
circ_datadf <- as.data.frame(circ_data)
circ_data$ctDNA.Dynamics <- NA #first we create the variable for the ctDNA & NAC combination, and we assign values
circ_data <- circ_data %>%
mutate(ctDNA.Dynamics = case_when(
ctDNA.Base == "POSITIVE" & ctDNA.postTx == "NEGATIVE" ~ 1,
ctDNA.Base == "POSITIVE" & ctDNA.postTx == "POSITIVE" ~ 2
))
circ_data$ctDNA.Dynamics <- factor(circ_data$ctDNA.Dynamics, levels=c("1","2"), labels = c("Clearance", "No Clearance"))
circ_data$Inclusion.status <- factor(circ_data$Inclusion.status, levels = c("Loco-regional", "Metastatic"))
circ_data$Location <- factor(circ_data$Location, levels = c("Oropharynx/Oral", "Larynx", "Unknown"))
circ_data$p16.status <- factor(circ_data$p16.status, levels = c("Negative", "Positive"))
circ_data$CPS.Scorev2 <- factor(circ_data$CPS.Scorev2, levels = c("≥1", "<1"))
surv_object <- Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event)
cox_fit <- coxph(surv_object ~ ctDNA.Dynamics + Age + Inclusion.status + Location + p16.status + CPS.Scorev2, data=circ_data)
ggforest(cox_fit, data = circ_data, main = "Multivariate Regression Model for PFS", refLabel = "Reference Group")
test.ph <- cox.zph(cox_fit)
#Univariate PFS cox regression for variables included in MVA
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Base!="",]
circ_data <- circ_data[circ_data$ctDNA.Base=="POSITIVE",]
circ_data$ctDNA.Dynamics <- NA
circ_data <- circ_data %>%
mutate(ctDNA.Dynamics = case_when(
ctDNA.Base == "POSITIVE" & ctDNA.postTx == "NEGATIVE" ~ 1,
ctDNA.Base == "POSITIVE" & ctDNA.postTx == "POSITIVE" ~ 2
))
surv_object <-Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event)
circ_data$ctDNA.Dynamics <- factor(circ_data$ctDNA.Dynamics, levels=c("1","2"), labels = c("Clearance", "No Clearance")) #univariate for ctDNA clearance post-treatment
cox_fit <- coxph(surv_object ~ ctDNA.Dynamics, data=circ_data)
summary(cox_fit)
Call:
coxph(formula = surv_object ~ ctDNA.Dynamics, data = circ_data)
n= 25, number of events= 20
coef exp(coef) se(coef) z Pr(>|z|)
ctDNA.DynamicsNo Clearance 2.085 8.048 0.767 2.719 0.00655 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95
ctDNA.DynamicsNo Clearance 8.048 0.1243 1.79 36.19
Concordance= 0.672 (se = 0.049 )
Likelihood ratio test= 11.61 on 1 df, p=7e-04
Wald test = 7.39 on 1 df, p=0.007
Score (logrank) test = 9.78 on 1 df, p=0.002
cox_fit_summary <- summary(cox_fit)
#Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)
[1] "HR = 8.05 (1.79-36.19); p = 0.007"
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Base!="",]
circ_data <- circ_data[circ_data$ctDNA.Base=="POSITIVE",]
surv_object <-Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event)
cox_fit <- coxph(surv_object ~ Age, data=circ_data) #univariate for age
summary(cox_fit)
Call:
coxph(formula = surv_object ~ Age, data = circ_data)
n= 25, number of events= 20
coef exp(coef) se(coef) z Pr(>|z|)
Age -0.02566 0.97467 0.02392 -1.073 0.283
exp(coef) exp(-coef) lower .95 upper .95
Age 0.9747 1.026 0.93 1.021
Concordance= 0.512 (se = 0.079 )
Likelihood ratio test= 1.12 on 1 df, p=0.3
Wald test = 1.15 on 1 df, p=0.3
Score (logrank) test = 1.16 on 1 df, p=0.3
cox_fit_summary <- summary(cox_fit)
#Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)
[1] "HR = 0.97 (0.93-1.02); p = 0.283"
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Base!="",]
circ_data <- circ_data[circ_data$ctDNA.Base=="POSITIVE",]
surv_object <-Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event)
circ_data$cStage <- factor(circ_data$cStage, levels = c("0-II", "III-IV")) #univariate for Stage
cox_fit <- coxph(surv_object ~ cStage, data=circ_data)
summary(cox_fit)
Call:
coxph(formula = surv_object ~ cStage, data = circ_data)
n= 25, number of events= 20
coef exp(coef) se(coef) z Pr(>|z|)
cStageIII-IV 0.9241 2.5196 0.6372 1.45 0.147
exp(coef) exp(-coef) lower .95 upper .95
cStageIII-IV 2.52 0.3969 0.7226 8.785
Concordance= 0.585 (se = 0.047 )
Likelihood ratio test= 2.55 on 1 df, p=0.1
Wald test = 2.1 on 1 df, p=0.1
Score (logrank) test = 2.24 on 1 df, p=0.1
cox_fit_summary <- summary(cox_fit)
#Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)
[1] "HR = 2.52 (0.72-8.78); p = 0.147"
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Base!="",]
circ_data <- circ_data[circ_data$ctDNA.Base=="POSITIVE",]
surv_object <-Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event)
circ_data$Inclusion.status <- factor(circ_data$Inclusion.status, levels = c("Loco-regional", "Metastatic")) #univariate for Stage/disease status
cox_fit <- coxph(surv_object ~ Inclusion.status, data=circ_data)
summary(cox_fit)
Call:
coxph(formula = surv_object ~ Inclusion.status, data = circ_data)
n= 25, number of events= 20
coef exp(coef) se(coef) z Pr(>|z|)
Inclusion.statusMetastatic -0.1110 0.8949 0.4582 -0.242 0.809
exp(coef) exp(-coef) lower .95 upper .95
Inclusion.statusMetastatic 0.8949 1.117 0.3646 2.197
Concordance= 0.556 (se = 0.056 )
Likelihood ratio test= 0.06 on 1 df, p=0.8
Wald test = 0.06 on 1 df, p=0.8
Score (logrank) test = 0.06 on 1 df, p=0.8
cox_fit_summary <- summary(cox_fit)
#Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)
[1] "HR = 0.89 (0.36-2.2); p = 0.809"
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Base!="",]
circ_data <- circ_data[circ_data$ctDNA.Base=="POSITIVE",]
surv_object <-Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event)
circ_data$p16.status <- factor(circ_data$p16.status, levels = c("Negative", "Positive")) #univariate for p16 status
cox_fit <- coxph(surv_object ~ p16.status, data=circ_data)
summary(cox_fit)
Call:
coxph(formula = surv_object ~ p16.status, data = circ_data)
n= 24, number of events= 20
(1 observation deleted due to missingness)
coef exp(coef) se(coef) z Pr(>|z|)
p16.statusPositive 0.2431 1.2752 0.6402 0.38 0.704
exp(coef) exp(-coef) lower .95 upper .95
p16.statusPositive 1.275 0.7842 0.3636 4.472
Concordance= 0.491 (se = 0.027 )
Likelihood ratio test= 0.14 on 1 df, p=0.7
Wald test = 0.14 on 1 df, p=0.7
Score (logrank) test = 0.14 on 1 df, p=0.7
cox_fit_summary <- summary(cox_fit)
#Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)
[1] "HR = 1.28 (0.36-4.47); p = 0.704"
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Base!="",]
circ_data <- circ_data[circ_data$ctDNA.Base=="POSITIVE",]
surv_object <-Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event)
circ_data$Location <- factor(circ_data$Location, levels = c("Oropharynx/Oral", "Larynx", "Unknown")) #univariate for disease location
cox_fit <- coxph(surv_object ~ p16.status, data=circ_data)
Warning: Loglik converged before variable 2 ; coefficient may be infinite.
summary(cox_fit)
Call:
coxph(formula = surv_object ~ p16.status, data = circ_data)
n= 25, number of events= 20
coef exp(coef) se(coef) z Pr(>|z|)
p16.statusPositive 2.431e-01 1.275e+00 6.402e-01 0.380 0.704
p16.statusUnknown -1.815e+01 1.309e-08 6.641e+03 -0.003 0.998
exp(coef) exp(-coef) lower .95 upper .95
p16.statusPositive 1.275e+00 7.842e-01 0.3636 4.472
p16.statusUnknown 1.309e-08 7.641e+07 0.0000 Inf
Concordance= 0.526 (se = 0.041 )
Likelihood ratio test= 3.49 on 2 df, p=0.2
Wald test = 0.14 on 2 df, p=0.9
Score (logrank) test = 1.92 on 2 df, p=0.4
cox_fit_summary <- summary(cox_fit)
#Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)
[1] "HR = -18.15 (0.78-76405588.46); p = 0.64"
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Base!="",]
circ_data <- circ_data[circ_data$ctDNA.Base=="POSITIVE",]
surv_object <-Surv(time = circ_data$PFS.months, event = circ_data$PFS.Event)
circ_data$CPS.Scorev2 <- factor(circ_data$CPS.Scorev2, levels = c("≥1", "<1")) #univariate for CPS score
cox_fit <- coxph(surv_object ~ CPS.Scorev2, data=circ_data)
summary(cox_fit)
Call:
coxph(formula = surv_object ~ CPS.Scorev2, data = circ_data)
n= 25, number of events= 20
coef exp(coef) se(coef) z Pr(>|z|)
CPS.Scorev2<1 0.3866 1.4720 0.4905 0.788 0.431
exp(coef) exp(-coef) lower .95 upper .95
CPS.Scorev2<1 1.472 0.6793 0.5628 3.85
Concordance= 0.564 (se = 0.062 )
Likelihood ratio test= 0.59 on 1 df, p=0.4
Wald test = 0.62 on 1 df, p=0.4
Score (logrank) test = 0.63 on 1 df, p=0.4
cox_fit_summary <- summary(cox_fit)
#Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)
[1] "HR = 1.47 (0.56-3.85); p = 0.431"
#Multivariate cox regression for OS
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Base!="",]
circ_data <- circ_data[circ_data$ctDNA.Base=="POSITIVE",]
circ_datadf <- as.data.frame(circ_data)
circ_data$ctDNA.Dynamics <- NA #first we create the variable for the ctDNA & NAC combination, and we assign values
circ_data <- circ_data %>%
mutate(ctDNA.Dynamics = case_when(
ctDNA.Base == "POSITIVE" & ctDNA.postTx == "NEGATIVE" ~ 1,
ctDNA.Base == "POSITIVE" & ctDNA.postTx == "POSITIVE" ~ 2
))
circ_data$ctDNA.Dynamics <- factor(circ_data$ctDNA.Dynamics, levels=c("1","2"), labels = c("Clearance", "No Clearance"))
circ_data$cStage <- factor(circ_data$cStage, levels = c("0-II", "III-IV"))
circ_data$CPS.Scorev2 <- factor(circ_data$CPS.Scorev2, levels = c("≥1", "<1"))
surv_object <- Surv(time = circ_data$OS.months, event = circ_data$OS.Event)
cox_fit <- coxph(surv_object ~ ctDNA.Dynamics + Age + cStage + CPS.Scorev2, data=circ_data)
ggforest(cox_fit, data = circ_data, main = "Multivariate Regression Model for OS", refLabel = "Reference Group")
test.ph <- cox.zph(cox_fit)
#Multivariate cox regression for OS v2
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("EORTC HNSCC ICI Clinical Data.csv")
circ_data <- circ_data[circ_data$ctDNA.available=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Base!="",]
circ_data <- circ_data[circ_data$ctDNA.Base=="POSITIVE",]
circ_datadf <- as.data.frame(circ_data)
circ_data$ctDNA.Dynamics <- NA #first we create the variable for the ctDNA & NAC combination, and we assign values
circ_data <- circ_data %>%
mutate(ctDNA.Dynamics = case_when(
ctDNA.Base == "POSITIVE" & ctDNA.postTx == "NEGATIVE" ~ 1,
ctDNA.Base == "POSITIVE" & ctDNA.postTx == "POSITIVE" ~ 2
))
circ_data$ctDNA.Dynamics <- factor(circ_data$ctDNA.Dynamics, levels=c("1","2"), labels = c("Clearance", "No Clearance"))
circ_data$Inclusion.status <- factor(circ_data$Inclusion.status, levels = c("Loco-regional", "Metastatic"))
circ_data$p16.status <- factor(circ_data$p16.status, levels = c("Negative", "Positive"))
circ_data$CPS.Scorev2 <- factor(circ_data$CPS.Scorev2, levels = c("≥1", "<1"))
surv_object <- Surv(time = circ_data$OS.months, event = circ_data$OS.Event)
cox_fit <- coxph(surv_object ~ ctDNA.Dynamics + Age + Inclusion.status + p16.status + CPS.Scorev2, data=circ_data)
ggforest(cox_fit, data = circ_data, main = "Multivariate Regression Model for OS", refLabel = "Reference Group")
test.ph <- cox.zph(cox_fit)